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Summary: 

Two accelerated imaginary-time evolution methods are proposed for the computation of solitary 
waves in arbitrary spatial dimensions. For the first method (with traditional power normalization), 
the convergence conditions as well as conditions for optimal accelerations are derived. In addition, 
it is shown that for nodeless solitary waves, this method converges if and only if the solitary wave 
is linearly stable. The second method is similar to the first method except that it uses a novel 
amplitude normalization. The performance of these methods is illustrated on various examples. 
It is found that while the first method is competitive with the Petviashvili method, the second 
method delivers much better performance than the first method and the Petviashvili method. 



1 Introduction 



In the study of nonlinear wave equations, solitary waves play an important role. In certain cases 
(such as in integrable equations), solitary waves can be calculated analytically. But in a majority of 
other cases, analytical expressions for solitary waves are not available. Important examples which 
have arisen recently in the study of physical systems include optical solitary waves in periodic media 
[1, 2, 3] and nonlinear matter waves in Bose-Einstein condensates (see, e.g., [4, 5, 6]). In such cases, 
one relies on numerical techniques to determine the shapes of solitary waves. Several types of 
numerical methods have been proposed and used for this purpose, such as the Newton's iteration 
method [7, 8], the shooting method [9], the nonlinear Rayleigh-Ritz iteration method [10], the 
Petviashvili method [11, 12, 13, 14], the imaginary-time evolution method [4, 5, 6, 15, 16, 17], and 
the squared-operator iteration methods [18]. The imaginary-time evolution method is attractive 
for its simple implementation, insensitivity to the number of dimensions, and high accuracy (due 
to its compatibility with the pseudo-spectral method). In addition, if it converges, it usually does 
so faster than the squared-operator iteration methods. 



The idea of the imaginary-time evolution method (ITEM) as applied to linear equations is quite old 
(see, e.g., [19, 20]). In the past decade, this method has also been applied to nonlinear equations 
[4, 5, 6, 15, 17]. In this method, one seeks the stationary solution of an evolution equation (usually, 
of the parabolic type) by numerically integrating that equation where time t is replaced by it (hence 
the name 'imaginary-time'), and normalizing the solution after each step of time integration to have 
a fixed L2 norm (called power by physicists). For linear equations, this method has long been known 
(see, e.g., [20]) to be equivalent to the problem of minimizing the energy functional of the physical 
system under the constraint that the solution being sought has a given power. Recently, this same 
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statement was shown to hold for nonlinear equations as well [17]. In [16], the authors treated the 
ITEM as a normalized gradient flow and proved its energy diminishing property. The ITEM, in its 
original form, is quite slow. In addition, it does not always converge to a stationary solution even 
if the initial function is quite close to the solution. In an effort to improve the convergence rate 
of the ITEM, the authors of Ref. [15] demonstrated that if the Sobolev gradients are used in the 
minimization of the energy functional, the convergence of the ITEM can be greatly accelerated (an 
equivalent possibility was mentioned in passing in [20], although no related details were provided 
there.) Alternatively, the authors of [17] used the steepest descent technique in the minimization 
of the energy functional and achieved fast convergence as well. However, these earlier studies or 
applications of the ITEM did not consider the conditions under which the ITEM and its accelerated 
versions would converge. So it was not clear when those methods could be used. In addition, the 
important practical question of establishing the conditions for the optimal acceleration of the ITEM 
was not considered either. 

An important question in the studies of solitary waves is their linear stability [21, 22, 23, 24, 25]. 
Most numerical techniques used to find solitary waves, such as the Newton's iteration method, 
the shooting method, and the Petviashvili method, yield no information about the stability of 
the solitary wave being obtained [7, 9, 12]. A remarkable fact about the ITEM and its properly 
accelerated version (with the usual power normalization), as we will show in this paper, is that the 
convergence of this numerical method is directly related to the linear stability of the corresponding 
solitary wave, provided that this wave is nodeless. This means that both existence and linear 
stability of these solitary waves can be obtained by this single numerical procedure. 

In this paper, we propose two new accelerated ITEMs for the computation of solitary waves in gen- 
eral nonlinear wave equations. Our acceleration technique is to introduce an acceleration operator 
to the imaginary-time equation, analogous to the preconditioning technique for solving systems of 
linear equations. For the first method, which uses power normalization, three important theoret- 
ical results are derived. One result is that convergence conditions of this method are explicitly 
obtained. This puts the application of this method on a solid theoretical footing. These conver- 
gence conditions show that in most cases, this method converges when the underlying solitary wave 
is nodeless, but there also exist cases when the method converges to solitary waves with nodes. 
Another result is that for nodeless solitary waves, this method converges if and only if the solitary 
wave is linearly stable. This connection between convergence of this method and linear stability 
of the underlying solitary wave is a novel property of this numerical method which has not been 
seen in other schemes. The third result is that explicit conditions for optimal acceleration of this 
method are obtained. These results provide the optimal practical implementation of this acceler- 
ated ITEM. The performance of this method is illustrated on various examples, and it is found to 
be competitive with the Petviashvili method. The second accelerated ITEM which we propose is 
similar to the first method except that it uses a non-traditional amplitude normalization. We will 
show through examples that this second method delivers better performance than the first method 
and the Petviashvili method. 

This paper is structured as follows. In Sec. 2, we introduce the original imaginary-time evolution 
method. In Sec. 3, we propose the first accelerated imaginary-time evolution method and derive its 
convergence conditions. In Sec. 4, we show that for nodeless solitary waves, the convergence of this 
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first method is directly linked to the linear stability of the solitary wave. In Sec. 5, we establish 
explicit conditions for the optimal acceleration of the first method. In Sec. 6, we propose the 
second accelerated imaginary-time evolution method, which employs the amplitude normalization. 
In Sec. 7, we apply both methods as well as the Petviashvili method to several examples, and show 
that the second accelerated imaginary-time evolution method delivers the best performance, while 
the first method is comparable to the Petviashvili method in performance. Sec. 8 concludes the 
paper. In the appendix, we attach a matlab code for one of the examples. 



2 Preliminaries on the original imaginary- time evolution method 



The problem we are interested in is the numerical determination of solitary waves in general scalar 
nonlinear wave equations in arbitrary spatial dimensions. To maintain the focus of the presentation, 
we first develop the theory for the iV-dimensional generalized nonlinear Schrddinger equation with 
an arbitrary potential. The extension of this theory to more general scalar nonlinear wave equations 
will be presented at the end of Sec. 5. 

The iV-dimensional generalized nonlinear Schrddinger equation with an arbitrary potential has the 
following form: 

iU t + V 2 U + F(\U\ 2 ,x)U = 0, (1) 
where x = (x±,X2, ■ ■ - xn) is a iV-dimensional spatial variable, 



V 2 - — + — > — (2) 

dx\ dx\ dx 2 N 



is the iV-dimensional Laplacian, and F(., .) is a real-valued function. This system is Hamiltonian. 
Solitary waves of Eq. (1) are sought in the form 

U(x,t) =tt(x)e i ' it , (3) 

where u(x) is a real-valued, localized function, and fi is a real parameter called the propagation 
constant. Then, from Eqs. (1) and (3), function u(x) is found to satisfy the equation 

L 00 u = fiu, (4) 

where 

L 00 = V 2 + F(n 2 ,x). (5) 



Equation (4) admits solitary waves for a large class of functions F(-u 2 ,x) [26, 27]. In this paper, 
we always assume that the solitary wave we are trying to obtain numerically does exist. 

In the original imaginary-time evolution method, one numerically integrates the equation 

u t = L 0Q u, (6) 
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which is obtained from Eq. (1) by replacing t with it (hence the name 'imaginary-time'), and then 
normalizes the solution after each step of time integration to have a fixed power. The power P of 
the solitary wave n(x) is defined as 

/CO 
u 2 (x;/x)eZx. (7) 
-co 

The simplest implementation of numerical time integration is to use the Euler method, whereby 
the ITEM scheme is: 



Un+l 

and 



P 



(u n+1 ,u n+1 ) 



u n +i, (8) 



u n+ i =u n + [L 00 u] u=Un At. (9) 
Here u n is the solution after the nth iteration, 

/CO 
/(x)*«7(x)dx (10) 
-co 

is the standard inner product in the iV-dimensional space of square-integrable functions, and the 
superscript "*" represents complex conjugation. Note that step (8) of the ITEM scheme guarantees 
that the power at every iteration is conserved: 

(u n ,u n ) = P, n = l,2, .... (11) 

Thus, if iterations (8)-(9) converge to a solitary wave it(x), then this u(x.) must satisfy Eq. (4), 
with its power being P and the propagation constant (i being equal to 

V = -js( L oo u ,u). (12) 

In the following sections, we will use two linear operators, whose explicit forms are given by the 
following definition. 

Definition 1 Operators Lq and L\ art defined as 

L = L 00 - \i = V 2 + F(u 2 , x) - n, (13) 

and 

Li = V 2 + F{u 2 , x) + 2u 2 F u 2 (u 2 , x) - n, (14) 

where F u 2 = dF/du 2 . 

Under these notations, Lqu = 0, and L\ is the linearization operator of Lqu with respect to u. 
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3 An accelerated imaginary- time evolution method 



A major drawback of the original ITEM (8)-(9) is that its convergence is quite slow, because the 
time step At has to be very small in order for it to converge (see Note 1 below for explanation). To 
overcome this difficulty, one idea is to use implicit time-stepping methods (such as the backward 
Euler method) to solve the imaginary-time equation (6) (see, e.g., [16]). Implicit methods allow 
larger time steps without causing divergence to the iterations. However, for nonlinear equations or 
in high dimensions, implicit schemes are difficult to implement, and their accuracies are often low 
(if finite-difference discretization is used). Another idea of accelerating the ITEM (see [15]) is to 
use the Sobolev gradients in the minimization of the energy functional. This idea is analogous to 
the preconditioning technique applied to the imaginary-time equation (6). Below we will extend 
this idea and propose a new accelerated ITEM that is explicit but fast-converging. 

In our accelerated ITEM, instead of evolving the original imaginary-time equation (6), we evolve 
the following "pre-conditioned" imaginary-time equation 

u t = M^ 1 [L 00 u - fiu] , (15) 

where M is a positive-definite and self-adjoint "preconditioning" operator. The stationary solution 
of this equation is still it(x). Applying the Euler method to this new equation, the new accelerated 
ITEM method (AITEM) we propose is: 



Un+l = Un+1\ -p. T , (16) 

\/ {u n+1 ,U n+1 } 

u n+1 =u n + M' 1 (L 00 u - ^u) u=Un At, (17) 



and 

(L 00 n, M-' l u) 



(u,M- 



(18) 



where P is the power defined in Eq. (7), which is pre-specified. Notice that our updating formula 
(18) for n n is quite special, different from the usual formula (12). This special updating formula 
(18) enables us to derive the convergence properties of the above AITEM, which we will do in this 
section. If M is the identity operator, then the above scheme is closely related to the original ITEM 
(8)-(9), and both have similar (slow) convergence properties. But if M takes other sensible forms 
[such as (45) below], convergence of the above AITEM will be much faster. In this paper, we will 
call M the acceleration operator. 



Before discussing the convergence conditions of the AITEM (16)- (18), we first present an example 
to demonstrate the drastic improvement in convergence of the AITEM (16)-(18) over that of the 
ITEM (8)-(9). 

Example 1 The nonlinear Schrodinger equation in one spatial dimension, 

u xx + u 3 = liu (19) 
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admits the solitary wave 

u(x) = v^sechx (20) 

with P = 4 and \x = 1. To apply the AITEM for this solution, we take M = 1 — c? xx , At = 1.5, and 
the Gaussian initial condition uq(x) = e~ x . The x-interval is taken to be [—15, 15], discretized by 
128 grid points, and the discrete Fourier transform is used to calculate u xx and to invert operator 
M. For the ITEM, the scheme parameters are the same, except that At = 0.01 (when At > 0.011, 
the ITEM diverges) . Defining the error as the L2 norm of the solution difference between successive 
iterations, i.e. [J(u n — u n _i) 2 dx] 1//2 , we find that for the error to drop below 10 -10 , the AITEM 
and ITEM take, respectively, 33 and 2160 iterations. Thus, the AITEM for this case is about two 
orders of magnitude faster than the original ITEM. The main reason for this drastic improvement 
of convergence is that for the AITEM, the time step At can be taken much larger (which is 1.5 
above) without causing divergence. This is made possible by the introduction of the acceleration 
operator M. 

Now we derive the convergence conditions for the AITEM (16)-(18). To do so, we introduce the 
following assumption on the kernel of L\ which holds in generic cases (in rare cases, this assumption 
may break down, see Fig. 7 in [28]). 

Assumption 1 If function F in Eq. (4) does not depend explicitly on certain spatial coordinates 
{xj 1 ,Xj 2 , ...,Xj k } (1 < ji, j'2, ■ ■■,jk < N), we assume that the only eigenf unctions in the kernel of 
L\ that are orthogonal to u(x) are the k translational-invariance modes du/dxj 3 (x),l < s < k. 

We also introduce a notation: for any operator L, we denote the number of its positive eigenvalues 
as p(L). 

Under the above assumption and notation, we have the following theorem on the convergence of 
the AITEM (16)-(18). 

Theorem 1 Let Assumption 1 be valid. Define At max = — 2/A m j n; where A m j„ is the minimum 
(negative) eigenvalue of operator C in Eq. (24). Then if At > At max , the AITEM (16) -(18) 
diverges. If At < At max , the following convergence statements on this AITEM hold. 

1. Ifp{L-i) = and P'(n) / 0, then the AITEM converges. 

2. Suppose p(Li) = 1, then p{M~ l L\) = 1. Denote the eigenfunction of this single positive 
eigenvalue of M~ X L\ as ip(x). Then if (ip,u) 7^ 0, the AITEM converges for P\n) > and 
diverges for P'(fJ-) < 0. If, however, (ip,u) = 0, the AITEM diverges. 

3. Ifp{L{) > 1, the AITEM diverges. 

Proof To analyze the convergence properties of the AITEM (16)- (18), we use the linearization 
technique. Let 

u n = u + u n , \u n \ < (21) 
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where u n {x) is the error. When this equation is substituted into the power-normalization step (16) 
and only terms of 0(u n ) retained, one obtains that the error is orthogonal to u(x): 

(u n ,u)=0, for all n. (22) 

Substituting Eq. (21) into (17)-(18) and linearizing, we find that the error satisfies the following 
iteration equation 

u n+l = (1 + At£)u n , (23) 

where operator C is 

y {u,M~ l u) J v ' 

Convergence of the AITEM depends on the eigenvalues of operator C. In view of the orthogonality 
constraint (22), the eigenvalue problem for £ that we need to consider is 

£^ = A* , * G S, (25) 

where 

5 = {*(x) : = 0}. (26) 

Note that when A ^ 0, by taking the inner product between equation = A\£ and u, one easily 
gets u) = 0, i.e. f £ S. Thus f G 5 in Eq. (25) constitutes a constraint only for zero 
eigenvalues of £. As we will show below, all eigenvalues A of £ are real, and all eigenfunctions of 
C form a complete set in S. Thus, in view of Assumption 1, the necessary and sufficient conditions 
for the AITEM to converge are that (i) for all non-zero eigenvalues A of £, 

-Kl + AAf<l; (27) 

(ii) for the zero eigenvalue of Eq. (25) (if exists), its eigenfunctions must be translational-invariance 
eigenmodes u Xj (which lead to a spatially translated solitary wave of Eq. (4) and do not affect the 
convergence of iterations). If At > At max with At max given in Theorem 1, then the left inequality 
in (27) is not met, thus the AITEM diverges. If At < At max , the left inequality in (27) is satisfied, 
hence we only need to consider the right inequality in (27), and the second condition (ii). 

We consider the condition (ii) first. Suppose Eq. (25) has a zero eigenvalue with eigenfunction 
\j>(x) j e 

{u,M- l u) v ' 

Differentiating Eq. (4) with respect to parameter ^, we find that 

LiUf, = u. (29) 

Hence the solution of Eq. (28) is au^, where a is a constant, plus functions in the kernel of L\. 
In view of Assumption 1 on the kernel of L\, we see that (^,u) = a(u^,u) = aP'(fi)/2. Thus if 
P'(lj) 7^ 0, then in order for $ £ S, a must be zero, hence ^ is in the kernel of L\. According 
to Assumption 1, ^ then must be a translational-invariance eigenmode which does not affect the 
convergence of iterations. 
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Next we consider the right inequality in (27), which is simply A < 0. We will analyze when this 
condition is met. To facilitate the analysis, we introduce the following new variables and operators 



* = M 1 / 2 ^, u = M- X l 2 u, Li = M- 1 / 2 L 1 M~ 1 /' 1 . (30) 
Then the eigenvalue problem (25) becomes the following equivalent one with eigenvalues unchanged: 

Ci> = A#, 4> £ S, (31) 

where 

t* = L 1 *-Hu, H = (Ll -1 5 , (32) 

(«,«) 

and 

5 = {#(x):(*,n)=0}. (33) 

It is easy to see that operator £ is self-adjoint in the set S due to L\ and M being self-adjoint and 
M being positive definite. Thus all eigenvalues A of £ are real, and all eigenfunctions ^ of £ form 
a complete set in S. As a result, all eigenvalues A of £ are real, and all eigenfunctions of £ form a 
complete set in S. In addition, since L\ is similar to M~ X L\ and in view of the Sylvester inertia 
law (see, e.g., Theorems 4.5.8 and 7.6.3 in [29]), p(L{) = p(M~ 1 L 1 ) = p(L 1 ). 

The eigenvalue problem (31) is equivalent to the one arising in the linear stability analysis of 
nodeless solitary waves (see, e.g., [2, 21] and, in particular, Eq. (2.3.9) in [2]), and has been well 
studied. To determine the sign of eigenvalue A, we expand ^ into the complete set of eigenfunctions 
of the self-adjoint operator L\ as 

*(x) = ]T 6 fc ^(x) + ( 6(A)^(x; X)dX, (34) 
k Jl 

where V>fc(x) an d V>(x; A) are the properly normalized discrete and continuous eigenfunctions of L\ 
with eigenvalues A& and A. The coefficients in Eq. (34) are given by bu = (i^k,^) and b = 
Function u can be expanded in a similar way with coefficients = u) and c(A) = (ip,u). When 
H = 0, the analysis is trivial and the convergence conditions in Theorem 1 can be easily obtained. 
Thus we only consider the H / case below. In this case, substituting expansions of u and ^ into 
Eq. (31), one finds that bu = Hck/{\k — A), 6(A) = Hc(X)/(X — A). Substituting these relations 
into the orthogonality condition ($f,u) = 0, we get 

\2 



If p(L\) = p{L\) = 0, i.e. all eigenvalues of L\ are negative, then Q(A) does not change sign when 
A > 0, hence equation (35) has no positive roots, and the right inequality in (27) is met. We 
have shown above that when P'(n) / 0, condition (ii) is met also, thus the AITEM converges. 
If p{L\) = p(L>i) > 1, denote Li's two positive eigenvalues as Ai and A2- If the corresponding 
expansion coefficients C\ and C2 in u are such that c\C2 ^ 0, then since Q(A) is continuous and 
monotonic in the interval (Ai, A2), and it approaches — 00 and 00 at the two ends of this interval, 
Q(A) obviously has a positive root between Ai and A2, hence the right inequality in (27) is not met, 
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and the AITEM diverges. If = (k = 1 or 2), then A = A& > is an eigenvalue of Eq. (31) with 
\P = ij)k, hence the AITEM also diverges. 

Now we consider the case p(L\) = p{L\) = 1. Denote this single positive eigenvalue of L\ as Ai, 
and its eigenfunction as if). If {ip, u) = 0, then A = Ai > is an eigenvalue of Eq. (31) with i$> = tp, 
hence the AITEM diverges. Below we consider the {ip,u) / case. Clearly Q(A) cannot have 
zeros when A > Ai. Whether Q(A) has a positive zero in (0, Ai) depends on the sign of Q(0)/H: if 
Q(0)/H < 0, then Q(\) has a positive zero, and vise versa. From Eqs. (31) and (32), we see that 
Q(0)/H = (L^ l ii,u). In view of Eq. (30), it follows that Q(0)/H = {L^u,u). Then from Eq. 
(29), we get Q(0)/H = P'{n)/2. Consequently, if P'(fi) > 0, the right inequality in (27) is met, 
and the AITEM converges; if P'{n) < 0, the AITEM diverges. Lastly, we notice from Eq. (30) 
that when ijj is an eigenfunction of Li, if) = M~ 1 / 2 ?/> is an eigenfunction of M~ 1 L± at the same 
eigenvalue, thus {tp,u) = (4>,u). This completes the proof of Theorem 1. □ 

Note 1 If M is taken as the identity operator, then our AITEM (16)-(18) becomes similar to the 
original ITEM (8)-(9). In this case, the smallest eigenvalue of C is A m j„ = — oo, which leads to 
At max = 0. However, in any computer implementation of this method, space is discretized. For the 
discretized operator of C, A m j n is finite, not — oo, thus At max > 0, hence this method (as well as 
the original ITEM) can still be used (see Example 1). But for any accurate spatial discretizations, 
h-min is large negative, hence At max is very small, which results in the slow convergence of this 
method and the original ITEM. For a more sensible choice of M such as (45) below, the smallest 
eigenvalues of C and its discretized version are both 0(1) and almost identical. This makes At max 
much larger, hence the AITEM (16)-(18) converges much faster (see Example 1). 

Theorem 1 is one of the main results of this article. It puts the application of the AITEM (16)-(18) 
on a firm theoretical basis. A corollary of this theorem can be readily established below. 



Corollary 1 Consider Eq. (4) with lim^^^ F(0, x) = 0. Under Assumption 1 and restriction 
At < At max , where At max is given in Theorem 1, the AITEM (16)- (18) diverges if the solitary wave 
it(x) has nodes (where u(x) = 0) and F u 2 (u 2 , x) > for ally:; z/u(x) is nodeless and F u 2 (n 2 , x) < 
for all x ; then the AITEM converges. 

To prove this corollary, we need the following lemma. 

Lemma 1 Consider the N -dimensional, linear Schrodinger eigenvalue problem 

V 2 ip - V(x)V> = \il>, (36) 

where V(x) — > as |x| — > oo. Let its discrete eigenvalues be arranged in the decreasing order: 
Ai > A2 > • • • > A m > 0, with the continuous spectrum being at X < 0. Let Vi(x) and V^(x) be 
two potentials such that V^(x) < Vi(x) for all x, and V^(x) ^ Vi(x). Then the discrete eigenvalues 
{A^} of Eq. (36) with potential V2 are larger than the corresponding eigenvalues {A^} of the 
same equation with potential V\, i.e., A^ > A^, k = 1,2, — 
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This lemma says that deepening the potential (making V(x) larger) shifts the eigenvalues of Eq. 
(36) downward. 



Proof of Lemma 1: Define the potential function 



y(x; a) = Vi(x) + a[V 2 (x) - Vi(x)], 



(37) 



where a is a real parameter. As a increases from to 1, V(x; a) changes from Vi(x) to V2(x). We 
now analyze how a discrete eigenvalue A of Eq. (36) changes as a continuously varies. For this 
purpose, we differentiate Eq. (36) with respect to a, and obtain 



In order for this equation to have a localized solution di^/da, its right-hand side must be orthogonal 
to the homogeneous solution ^(x; a). This yields the relation 



According to our assumption, V^(x) < Vi(x) for all x, and V^(x) / Vi(x); thus dX/da > 0. This 
means that A^ > X^\ k = 1, 2, — Hence Lemma 1 is proved. □ 

With this lemma, we now prove Corollary 1. 

Proof of Corollary 1: First, we recall the definitions of the two operators Lq and L\ in Eqs. 
(13) and (14), and the fact of Loti(x) = 0, i.e., Lq has a zero eigenvalue A a = with eigenfunction 
u(x). We will also use a well known result about linear Schrodinger operators, which says that the 
largest eigenvalue of those operators is simple and the corresponding eigenfunction is nodeless [30]. 
We now consider two possibilities in regards to the nodes of u(x). 

(1) Suppose u(x) has at least one node. Then, by the aforementioned property of linear Schrodinger 
operators, Lq must also have at least one positive eigenvalue Aft > whose eigenfunction is nodeless. 
Comparing the two Schrodinger operators Lq and L\, we see that the difference in their potentials 
is 2w 2 F u 2(u 2 ,x). If F u 2(n 2 ,x) > and F n 2(-u 2 ,x) ^ 0, then it is seen from Lemma 1 that operator 
L\ must have at least two positive eigenvalues corresponding to the eigenvalues A a and A& of Lq. 
Hence according to Theorem 1, the AITEM diverges. 

(2) If n(x) is nodeless, then zero is the largest eigenvalue of Lq. When F n 2(-u 2 ,x) < and 
F u 2(n 2 ,x) ^ 0, then by Lemma 1, the eigenvalues of L\ are all negative. Thus according to 
Theorem 1, the AITEM converges. This completes the proof of Corollary 1. □ 

Corollary 1 can be readily used on certain equations, and two examples are shown below. 

1. Consider the focusing nonlinear Schrodinger equation with a single- well potential: 




(38) 



dX _ ((V 1 -V 2 )M) 
da {ip, ip) 



(39) 



u X x + 6sech 2 x u + u 3 = fiu. 



(40) 
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This equation admits a family of solitary waves with nodes [see an example in Fig. 1(a)]. Here 
F u 2(u 2 ,x) = 1 > 0, thus according to Corollary 1, the AITEM diverges for these solutions, which 
we have confirmed numerically. 



2. Consider the defocusing nonlinear Schrodinger equation with a single-well potential: 



u xx + 6sech 2 x u — u 3 = 



fiU. 



(41) 



This equation admits a family of nodeless solitary waves [see an example in Fig. 1(b)]. Here 
F u 2(u 2 ,x) = — 1 < 0, thus according to Corollary 1, the AITEM converges for these solutions 
under restriction At < At max (even though it can be verified that P'(n) < here). 

Corollary 1 suggests that the AITEM tends to diverge for solitary waves with nodes (such solutions 
are often called "excited states"). Indeed, this if often the case. However, there are examples where 
the AITEM converges for solitary waves with nodes. This could occur if all the eigenvalues of L\ are 
negative. For a solitary wave u(x) with nodes, since Lqu = 0, Lq must also have positive eigenvalues 
[30]. But if F u 2 (u 2 , x) < for all x, then the eigenvalues of L\ are lower than those of Lq according 
to Lemma 1. If the potential difference \2u 2 F u 2 (u 2 , x)| between Lq and L\ is large enough, it is 
possible for all the eigenvalues of L\ to be pushed below zero, resulting in the convergence of the 
AITEM. One such example is given below. 

Example 2 Consider the defocusing nonlinear Schrodinger equation with a double-well potential: 



For this equation, F u 2(u 2 ,x) = — 1 < 0. This equation admits a family of single- node solitary 
waves, whose power curve is displayed in Fig. 2(a). It is seen that P'(fi) < for the entire family. 
Two representative solutions with powers P = 3 and P = 10 are displayed in Fig. 2(b). The 
spectra of operator L\ for these two waves are plotted in Fig. 2(c, d), respectively. At the lower 
power P = 3, p{L\) = 1. However, at the higher power P = 10, \2u 2 F u 2 (u 2 , x)| = 2-u 2 is large 
enough such that p(L\) = 0, hence according to Theorem 1, the AITEM converges. 



4 Connection between convergence and linear stability 

The convergence theorem 1 strongly resembles the linear stability conditions of nodeless solitary 
waves in the generalized NLS equations (1) [2, 21, 23, 24, 31]. For such solitary waves, the following 
stability conditions have been established [2]: 

For a nodeless solitary wave in Eq. (1), (i) ifp(L\) < 0, the wave is linearly stable; (ii) ifp(Li) = 1, 
the wave is linearly stable if P'(fi) > and unstable if P'((i) < 0; (Hi) if p(L{) > 1, the wave is 
linearly unstable. 

This stability result is almost identical to our convergence theorem 1, indicating that the AITEM 
(16)-(18) usually converges to a nodeless solitary wave if and only if this wave is linearly stable. 



Uxx + 6 sech 2 (x + 1) + sech 2 (x — 1) u — v? = fiu. 



(42) 
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The only notable difference between the above stability and convergence results is in the case 
of p(Li) = 1, where the convergence theorem has a condition on (ij),u) which is absent in the 
stability theorem. However, for nodeless solitary waves u(x) with p(L{) = 1, condition (if),u) / 
in Theorem 1 is met in generic cases, thus the stability and convergence results are the same for 
this case as well. For the following two choices of the acceleration operator M, we can actually 
show that condition (ij),u) ^ is strictly satisfied. 

The first choice is when M is the identity operator, where the AITEM becomes similar to the 
original ITEM (8)- (9). In this case, if p{L\) = 1, then VK X ) is the eigenfunction of the largest 
eigenvalue of the Schrodinger operator Li, which is known to be nodeless [30]. Since u(x) is 
nodeless as well, (tp, u) thus is non-zero. 

The second choice is a very practical one, M = \i — V 2 , which will be shown to yield optimal 
acceleration for a large class of equations in the next section. For this M, the eigenvalue equation 
M~ X L\^) = Xtp can be rewritten as the following Schrodinger equation 

VV-^ + ^V' = 0, (43) 

where 

V(x) = F(u 2 ,x) + 2u 2 F u2 (u 2 ,x). (44) 

When p{L\) = p(M~ l L\) = 1 and A is the single positive eigenvalue of M~ l L\ [i.e. Eq. (43)], 
it is easy to show using the spectral properties of Schrodinger operators that the corresponding 
eigenfunction tp(x.) is nodeless (a similar fact for the ID case can be found in [32]). Since u(x) is 
also nodeless, (ip, u) / 0. 

It is remarkable that the AITEM (16)-(18) can not only produce solitary waves, but also determine 
their linear stability properties. This is like "killing two birds with one stone". To the authors' 
knowledge, there are no other numerical methods for solitary waves which possess this same prop- 
erty. Note, however, that this property holds only for nodeless solitary waves. For solitary waves 
with nodes, this connection between convergence of the AITEM and linear stability of the solitary 
wave can break down. 



5 Optimal acceleration of the imaginary- time evolution method 

In the AITEM (16)-(18) for Eq. (4), a practical choice of the acceleration operator M is 

M = c - V 2 , c> 0. (45) 

The reason for this choice is two fold: first, M~ 1 is very easy to compute by the fast Fourier 
transform; second, all eigenvalues of C are 0(1), which makes At max = 0(1) as well. For this M, 
an important question then is: at what value of c does the AITEM converge the fastest? We will 
answer this question in this section. Note that the use of a seemingly more general form of the 



12 



acceleration operator M = c± — C2V 2 is equivalent to (45) by a rescaling of the time step At in Eq. 
(17), and hence does not warrant consideration. 



First, we quantify the convergence rate of the AITEM. From Eq. (23), we see that the error u n 
evolves in proportion to R n , where 

R(c,At) = max{|l + A min At|, \l + A max At\}, (46) 

and A min (c), A max (c) are the smallest and largest non-zero eigenvalues of Eq. (25). If R < 1, the 
AITEM converges, and vise versa. In this section, we assume that the AITEM converges (under 
the stepsize restriction in Theorem 1). This means that both A m j n (c) and A max (c) are negative. 
The parameter R characterizes the rate of convergence and will be called the convergence factor. 
Smaller R leads to faster convergence. For fixed c, it is seen from Eq. (46) that the smallest R is 
reached at 

At = AU(c) = — , (47) 

whence 

r> / \ n/ \ 4. \ A m i n — A max Sao\ 

i?*(c) = R(c, At,) = — . (48) 



The value of c which makes R*(c) minimal gives optimal acceleration of the AITEM and will be 
denoted as c opt . The determination of c op t is the focus of this section. 

Before analytically determining c op t, we first present a numerical example. Consider the NLS 
equation (19) with \x = 1 again. For each c value in Eq. (45), we have numerically obtained 
Amin(c) and A max (c) of operator C by discretizing Eq. (25) and turning it into a matrix eigenvalue 
problem. The resulting i?*(c) function is then obtained from Eq. (48) and plotted in Fig. 3(a). 
We see that the minimum of R*(c) occurs at c = 1, thus c op t = 1. At this c value, dependence of 
the convergence factor R on the timestep At can be calculated from Eq. (46) and is displayed in 
Fig. 3(b). We see that when At > 2, R > 1, thus iterations diverge. When < At < 2, iterations 
converge, and fastest convergence occurs when At rs 1.51, which is the value from Eq. (47). 

In the above numerical example, it is observed that c op t = fi. Is this a coincidence? The answer 
is negative. Below, we will show that for a large class of equations (4) with localized potentials, 
c pt = This result is stated in the following theorem. 

Theorem 2 Consider Eq. (4) with Hm^^cc F(0, x) = 0. If V(x) given in (44) does not change 
sign, then c op t = (J, in the AITEM (16)- (18). If V(x) changes sign, then c = fj, is not optimal in 
the generic case. 

The generic case will be defined later in Lemma 3. It is noted that when lim^^^ F(0,x) = 0, 
H > for solitary waves in Eq. (4). 

To facilitate the proof of Theorem 2, we first establish a few lemmas. 
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Lemma 2 If lim^^^ F(0,x) = in Eq. (4), then the continuous spectrum of operator C in (25) 
is given by 

AeKil, for c >f i; 

(49) 

A G [-£,-1), for c< ii. 
Moreover, if (A,*) is a discrete eigenmode of Eq. (25), then 

IdA 1 , N 

< " a * = WJik < 7 < 50 > 



Proof: First, since C — > M _1 Li as x — ► oo, the continuous spectrum of £ is then the same as 
that of M~ 1 Li, which can be easily shown to be (49). Next, we consider how a discrete eigenvalue 
A changes with c. The eigenvalue equation = A* can be rewritten as 

Li* - Hu = AMf, (51) 

where H = (Li\l/,M -1 u)/(«, M -1 ti). Differentiating Eq. (51) with respect to c, then taking its 
inner product with \I> and noticing $ £ 5, we get 

<(Li-AM)^,tf) = <^Mtt + Atf,tt). (52) 

Since both Li and M are self-adjoint and utilizing Eq. (51), we have 

((Ll _ AM)—, vl/) = ( — , (Lx - AM)*) = H(—,u) = H-{% «>, (53) 



which is zero since $ 6 S. Thus from Eq. (52), we get 

IdA (*,*) 



A dc (*,M*)' 
Since M*) > c(*, *) > 0, Lemma 2 is proved. □ 



(54) 



Lemma 3 Consider Eq. (4) with ^m\ x \^oo ^(0, x) = 0. If either A m j n (/x) = —1 or A max (p,) = —1, 
then c op t = fJ" If neither A max (/j,) nor A m j„(/z) equals —1, then in the generic case where 



1 dA, n 



A r , 



dc 



1 dA r 



C=fl 



A r , 



dc 



(55) 



C=/l 



c = (i is not optimal. 



Proof: Differentiating formula (48) of i?*(c) with respect to c, one gets: 



dR* 



*2A m i n A max 



dc (A m j n + A max )^ 



1 dA r 



1 dA r 



A, 



dc 



dc 



(56) 
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The factor in front of the square brackets above is positive since both A max and A min are negative 
(see beginning of this section). Following the assumption of this lemma, suppose, for definiteness, 
that A min {/j,) = — 1. When c decreases from //, the lower edge of the continuous spectrum decreases 
as — fi/c [see Eq. (49)]. Due to inequality (50), all discrete eigenvalues of C decrease slower than 
—fj>/c, thus A mira (c) = — fi/c. Then using Eqs. (50) and (56) one obtains: 



dRif 

dc 



< 



m%n li -max 



{A min + 1 ^max I 



H c 



0, c < ji. 



(57) 



When c increases from fi, the lower edge of the continuous spectrum is always at —1 [see Eq. (49)], 
while discrete eigenvalues of C all increase [see Eq. (50)] , thus A m j n (c) = —1. Then using Eqs. 
(50) and (56) we get 



dR* 
dc 



> 



2A.7717" 71, A. 



min 1 *--ma:r 



[0 + 0] = 0, c> 11. 



{.A m in ~\~ A m ax . 

Inequalities (57)-(58) mean that R*{c) has a global minimum at c = fi, thus c op t 
— 1, following similar arguments one can show that c op t = /t as well. 



(i. 



If A, 



(58) 



If, however, A m j n (/x) < — 1 and A ma x{n) > — 1, then under the condition (55) which holds in the 
generic case, Ri{fJ>) exists and is not equal to zero. Thus the minimum of i?*(c) is not at c = fi, i.e. 
c = fi is not optimal. Lemma 3 is thus proved. □ 



In the following two lemmas, we establish the conditions under which either A m j n (//) 



-1 or 



A m ax{n) = — 1- For convenience, we define M$ = M c=a 
A m ax{tj) are the smallest and largest eigenvalues of Cq. 



/x - V 2 , £ = £ c =u- Then A min (/x) and 



Lemma 4 Suppose limix^^ F{0, x) = in Eq. (4), then Cq and MJ~ L\ both do not have con- 
tinuous spectrum, and their discrete eigenvalues accumulate at —1. In addition, if the smallest 
eigenvalue \ m in of M^ l L\ is —1, then A m j n (/i) = —1; if M^ l L\ has two or more eigenvalues that 
are less than —1, then A m j„( / u) < —1. Similarly, if the largest eigenvalue X m ax of Mq L\ is —1, 
then Amax{n) = —1/ if M$ 1 Li has two or more eigenvalues greater than —I, then A max {n) > —1. 



Proof: Operators Mq Li and Co do not have continuous eigenvalues, as follows from Eq. (49). 
The eigenvalue equation M Q ~ 1 Liip = At/> is the same as the Schrodinger equation (43). Using 
well known spectral properties of the Schrodinger operators, we know that operator M ~ 1 Li has an 
infinite number of (discrete) eigenvalues, which accumulate in such a way that (1 + A) -1 approaches 
either +00 or —00, or both (this fact for the ID case can be found in [32]). Thus the accumulation 
point of eigenvalues for Mq 1 L\ is —1. 

Applying the technique in the proof of Theorem 1 (see Eqs. (34)-(35), or [2, 21]) on the eigenvalue 
problem for £q, one can show that between every two adjacent eigenvalues of Mq 1 L\, there is an 
eigenvalue of Cq. Thus, —1 is also an accumulation point of eigenvalues for Cq. By the same reason, 
if the smallest eigenvalue \ m in for M ~ 1 Li is —1, then A m j n (^) = —1; if M ~ 1 L\ has two or more 
eigenvalues that are less than —1, then A m j n ( / u) < — 1. Similarly, if the largest eigenvalue X m ax f° r 
Mq 1 L\ is —1, then A max (/x) = —1; if Mq X L\ has two or more eigenvalues greater than —1, then 
Amax { lA > — 1- Lemma 4 is thus proved. □ 
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Lemma 5 Consider Eq. (4) with lim^^^ F(0, x) = 0. For operator Mq L 1; if V(x) > for all 
x, then its smallest eigenvalue \ m in is — 1; if V(x) < /or aZZ x ; then its largest eigenvalue \ m ax is 
— 1; i/V(x) changes sign, then there is an infinite number of its discrete eigenvalues on both sides 
of -I. 



Proof: Consider the eigenvalue equation M 1 Liip = which is the same as Eq. (43). Taking 
the inner product of Eq. (43) with ip, we get 

1 + A= (Mo^W- (59) 

If V(x) > for all x, then since M$ is a positive definite self-adjoint operator, the right hand 
side of the above equation is non-negative, thus A > —1. Due to Lemma 4, eigenvalues of Mq 1 L\ 
accumulate at —1, thus \ m in = — 1- By similar arguments, if V(x) < for all x, then X ma x = — 1- 
If V(x) changes sign, then Eq. (43) has infinite numbers of discrete eigenvalues on both sides of 
— 1 [32]. Hence Lemma 5 is proved. □ 



With these lemmas, we are now ready to prove Theorem 2. 



Proof of Theorem 2: If V(x) does not change sign, then by Lemmas 4 and 5, A m j n (/x) = — 1 
or A max (n) = —1, hence by Lemma 3, c op t = n- On the other hand, if V(x) changes sign, Lemma 
5 indicates that there are infinitely many eigenvalues of Mq 1 L\ both below and above —1. Then 
by Lemma 4, there are also infinitely many eigenvalues of Cq both below and above —1, so that 
A m in(^) < — 1 and A max (fi) > — 1. Then in the generic case defined by equation (55), c = fx is not 
optimal by Lemma 3. Theorem 2 is thus proved. □ 

In practical implementations of the AITEM, after c op t = n in Eq. (45) is chosen, one still needs to 
select the time step At. The best choice of At which leads to fastest convergence is given by Eq. 
(47). Since the exact values of A m j„ and A max are usually not available, below we give the interval 
of values where the optimal time step can be found. When V(x) > for all x, A m j n (/x) = — 1, 
and — 1 < A max (fi) < 0, hence 1 < At*(/x) < 2. When V(x) < for all x, A m j„(//) < — 1 and 
Am ax (ji) = -1, hence At*(/x) < 1. 

In some physical problems, the assumption limix^^ F(0, x) = is not met, i.e. the potential in 
Eq. (4) is not localized. One example is the following type of equations 

V 2 u + F(x)n + G(n 2 ,x)n = fiu, (60) 

where V(x) is a periodic function in x, and G(0,x) ->0asx^ oo. For these equations, if we take 
the acceleration operator in the form of M = c — V 2 — V(x), then we can also show that under 
conditions analogous to those in Theorem 2, c op t = \i. However, for this form of M, it is not easy 
to compute M _1 . So in practice, it may still be better to use the simple form (45) for Eq. (60) 
instead. In that case, c op t is not known analytically and may need to be estimated by trial and 
error (see Example 4 in Sec. 7). Our experience shows that in many cases, taking a suboptimal c 
does not slow down the method significantly as long as the c taken is not far away from c op t- Thus 
using the form (45) of M for Eq. (60) does not constitute a significant disadvantage. 
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Lastly, we would like to point out that most of the results in Sees. 3, 4 and 5 can be extended to 
a wider class of equations 

Vu + F(u 2 ,x)u = frn, (61) 



where u is a real-valued localized function, F is also real-valued, \i is a real parameter, and V is 
a general linear self-adjoint semi-negative-definite constant-coefficient pseudo-differential operator. 
The previous equation (4) is a special case of (61) with T> = V 2 . For this general equation (61), the 
AITEM is still (16)-(18), with V 2 replaced by V. The convergence conditions of this AITEM are 
still the same as those in Theorem 1, except that V 2 in the definition (14) of L\ is replaced by V. 
Regarding optimal acceleration, if M = c — T> is taken, then we can also show that c op t = n when 
limixi^oo F(0, x) = and V(x) as defined in (44) does not change sign. The connection between 
convergence of the AITEM and linear stability of a nodeless solitary wave can also be extended 
to other types of nonlinear wave equations. For instance, consider the generalized one-dimensional 
Korteweg de Vries (KdV) equation 



where U(x,t) and F(-) are both real- valued functions, F(0) = 0, and T> is a general linear self- 
adjoint semi-negative-definite constant-coefficient pseudo-differential operator. Looking for moving 
solitary waves of the form U (x, t) = u(x — fit), where jj, is a real parameter, we get an equation for 
u which is a special form of (61). For Eq. (62), it has been shown that if u(x) has no nodes and 
p{L\) = 1, then the solitary wave u(x — fit) is linearly stable if P'{(j) > and linearly unstable if 
P'(n) < [33, 34]. In this case, for generic acceleration operators M where the eigenfunction of 
M _1 Li's positive eigenvalue is non-orthogonal to u(x), the convergence conditions of the AITEM in 
Theorem 1 are the same as the stability conditions above. Thus the connection between convergence 
of the AITEM and linear stability of the solitary wave holds for Eq. (62) as well. 

6 The accelerated imaginary- time evolution method with ampli- 
tude normalization 

The AITEM (16)-(18) discussed above employed the power normalization [see (16)], which has 
commonly been used in all previous ITEM-type methods. A consequence of this normalization is 
that when P'(fi) = 0, this method is doomed to fail. In certain cases when P'(fi) < (see case 2 in 
Theorem 1), this AITEM diverges as well. In this section, we will propose a different normalization 
for the AITEM which can overcome the above difficulties. 

This new normalization is the amplitude normalization. In other words, instead of fixing the power 
of the solution we are seeking, we fix the amplitude of u(x) (i.e., the largest value of |u(x)|). Thus, 
this new AITEM we propose for Eq. (4) is 



U t + VU + F(U 2 )U =0, 



(62) 



J X 



A 



(63) 



Un+l = 



Un+l \ max 



U n+ l =U n + M 1 (L 00 U - flu) 



u=u„, [M=n„ 



At, 



(64) 
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and 



{u,M- l u) 



(65) 



where A = \u\ max is the pre-specified amplitude of u(x). This AITEM with amplitude normaliza- 
tion, which we denote as AITEM(A.N.), can converge regardless of the value of P'{y). 

We have tested this AITEM(A.N.) on various examples, and found that it is almost always more 
superior than the first AITEM (16)- (18). This superiority is reflected on two aspects: (i) in cases 
where the first AITEM does not converge, the AITEM (A.N.) often can converge; (ii) in cases 
where the first AITEM converges, the AITEM(A.N.) often converges faster. 

To illustrate the first aspect of the AITEM (A.N.) superiority, we consider the following example. 



Example 3. The 2D NLS equation 

U X X + Uyy + U 3 = JJLU (66) 

admits a family of single-hump (fundamental) solitary waves, all of which have the same power 
P = 11.70. Since P'{n) = everywhere, the first AITEM (16)- (18) clearly can not converge 
to these solitary waves. However, the AITEM (A.N.) can not only converge, but also converge 
very fast. To illustrate, we select the solitary wave with amplitude A = 1, whose corresponding 
propagation constant is \i = 0.2054. This solitary wave is displayed in Fig. 4(a). We take the 
acceleration operator M as (45), the initial guess as a Gaussian hump uo(x,y) = e~( x +y ^ 10 , the 
spatial domain as a square with side length of 30, with each dimension discretized by 128 points. 
We also take (c, At) = (0.5,1), which is nearly optimal. With these choices, the AITEM (A.N.) 
rapidly converges to the solitary wave [see the error diagram in Fig. 4(b)]. Indeed, it takes only 27 
iterations for the error to fall below 10~ 10 . For comparison, we applied the standard Petviashvili 
method to this example [11], taking the same spatial discretizations and initial condition as above. 
The error diagram of the Petviashvili method is also displayed in Fig. 4(b). We see that the 
AITEM (A.N.) is faster than the Petviashvili method by about 40%. 



The second aspect of the AITEM (A.N.) superiority over the first AITEM will be illustrated by 
examples in the next section. 

It should be pointed out that for the AITEM (A.N.), the connection between convergence of the 
scheme and linear stability of the solitary wave disappears. Regarding its convergence conditions, 
this question can be analyzed by techniques similar to those used in Sees. 3 and 5. This will not 
be done in this article, and will be left for future studies. 



7 Examples illustrating convergence rates of the AITEMs 



In this section, we will apply the two proposed AITEMs (16)-(18) and (63)-(65) to two physi- 
cal examples, and compare their convergence speeds. In addition, we will compare them to the 
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Petviashvili method whenever applicable. 



Example 4 Let us first consider a 2D NLS equation with a periodic potential, 

Uxx + u yy + Vo ^cos 2 x + cos 2 yj u + u 3 = fiu, (67) 

which has recently attracted much interest due to its application to optical lattices and Bose- 
Einstein condensates [3, 35, 36]. This equation admits a family of nodeless solitary waves. One 
of them with V = 3, P = 3, fi = 3.7045 and amplitude A = 1.0384 is displayed in Fig. 5 
(a). For this wave, p{L\) = 1 and P'(fi) > 0, thus the AITEM converges for generic choices of 
M. We applied the AITEM and AITEM(A.N.) to search for this solitary wave. The acceleration 
operator M was taken as (45), the spatial domain taken to be a square with side length of 107T, 
with each dimension discretized by 128 points. As the initial guess, we took a Gaussian profile 
uo(x,y) = exp(— x 2 — y 2 ). For our choice of (45), c op t is not known analytically. Hence we scanned 
values of c and At in the AITEMs and found that taking c = 0.7 and At = 1 yielded the fastest 
(or nearly fastest) convergence for both AITEM and AITEM(A.N.). The error diagrams versus the 
iteration number for these two methods are shown in Fig. 5(b). We see that an error below 10 -10 
is reached by the AITEM and AITEM (A.N.) in about 310 and 130 iterations respectively. Thus 
the AITEM (A.N.) converges much faster than the AITEM. In the appendix, the matlab code of 
AITEM(A.N.) on this example is attached, so that the reader can test this method themselves. 
For Eq. (67), the original Petviashvili method does not apply [11]. Several generalizations of that 
method for equations such as (67) have been proposed recently [13, 14, 37]. A comparison between 
those generalized Petviashvili methods and the AITEMs in this paper will be considered elsewhere. 

Example 5 Another example we consider is an equation whose linear part is not of the Schrodinger 
type. Specifically, we consider the integrable Kadomtsev-Petviashvili (KP) equation, whose sta- 
tionary solutions satisfy the equation 

Uxx — dx 2 Uyy + u n = fiu, n = 2 (68) 

and the constraint 



u(x, y)dy = . (69) 
The analytical expression for solitary waves of Eq. (68) is [38] 

3 + fi 2 y 2 — fix 2 



u — 12(i — 

(3 + fi 2 y 2 + fix 2 ) 2 ' 



(70) 



and their power function is P{fi) = 24^^///. It is known that p(L\) = 1 for solutions (70) [12], 
and P'(fi) > in view of the above power formula. Thus the AITEM converges to these localized 
solutions for generic choices of M (see end of Sec. 5). We applied the AITEM, AITEM (A.N.) 
and the Petviashvili method [11] to compute one of these solutions with fi = 1, whose profile is 
shown in Fig. 6(a). The spatial domain was taken as a square with side length of 120, discretized 
by 512 points in each dimension. The acceleration operator M in the AITEMs was taken as 
M = c — d xx + d~ 2 d yy . Since V(x,y) = 2u changes sign here, generically c op t ^ fi in the AITEM. 
For the AITEM and AITEM (A.N.), we took (c, At) = (1.4,1.7) and (1.4,1.4) respectively, which 
yield near-optimal convergence for the underlying methods. For all three methods, we imposed 
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the constraint (69) by setting the coefficient of the = Fourier harmonic to be zero at every 
iteration. The error diagrams for these methods are displayed in Fig. 6(b). This time, an error 
below 1(T 10 is reached by the AITEM, AITEM(A.N.) and the Petviashvili method in about 140, 40 
and 90 iterations respectively. Thus the AITEM (A.N.) converges much faster than the Petviashvili 
method, while the Petviashvili method converges faster than the AITEM. 

We have studied solitary waves in the generalized KP equation (68) (with n / 2) as well. In those 
equations, if n > 7/3, solitary waves are unstable [39]. In such cases, we found that the AITEM 
(16)-(18) diverged. Thus the connection between convergence of that AITEM and the stability of 
the underlying solitary wave holds for the generalized KP equation as well. 

From the above examples, we conclude that the AITEM (A.N.) converges faster than the AITEM 
and Petviashvili methods, while the AITEM and the Petviashvili method are comparable in con- 
vergence speeds. 

Before concluding the paper, we make a comment on the practical implementations of the AITEMs 
(and other methods such as the Petviashvili method as well). If the convergence theorem predicts 
that a method diverges for a solitary wave, sometimes iterations can still converge (at least to 
a certain accuracy). This can happen, for instance, when the iteration operator C has a single 
unstable symmetric eigenmode, but the initial condition is chosen to be, and the final solution is, 
strictly anti-symmetric. In that case, the unstable symmetric eigenmode may not yet be excited 
before the iterations have already converged. Such an example is Eq. (42), which admits anti- 
symmetric solitary waves. When P = 3 (see Fig. 2), p(L\) = 1, and P'(S) < 0. Hence the AITEM 
(16)-(18) should diverge. However, we found that if we used anti-symmetric initial conditions, then 
the AITEM iterations can converge to the solution with error below 10~ 12 . However, if a small 
symmetric component was introduced into the initial condition, the iterations would diverge. Thus, 
the convergence results obtained in this paper must be understood as pertaining to generic initial 
conditions. If non-generic initial conditions are used in practical implementations of the AITEMs, 
better convergence behavior may be observed. 



8 Summary 



In this paper, we proposed two accelerated imaginary-time evolution methods for computations 
of solitary waves in arbitrary spatial dimensions. The first method is the AITEM (16)-(18) with 
the conventional power normalization. For this method, the convergence conditions were derived. 
These conditions show that this method usually converges to nodeless solutions, but there also exist 
cases when this AITEM converges to solitary waves with nodes. For nodeless solutions, we also 
showed that this AITEM converges if and only if the solitary wave is linearly stable. Conditions 
for optimal accelerations of this AITEM were also derived. The second method we proposed is 
the AITEM (A.N.) (63)-(65) which employs a novel amplitude normalization. Both methods were 
applied to several examples of physical interest, and we found that the AITEM (A.N.) delivers the 
best performance, while the AITEM and the Petviashvili method are comparable in performance. 
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Appendix: Matlab code of AITEM (A.N.) for Example 4 

Lx=10*pi; Ly=10*pi; N=128; c=0.7; DT=1; 
max_iteration=10000; error_tolerance=le-10; 

dx=Lx/N; x=-Lx/2:dx:Lx/2-dx; kx=[0: N/2-1 -N/2 : -1] *2*pi/Lx; 
dy=Ly/N; y=-Ly/2:dy:Ly/2-dy; ky=[0: N/2-1 -N/2 : -1] *2*pi/Ly ; 
[X,Y]=meshgrid(x, y) ; [KX , KY] =meshgr id (kx , ky) ; 

A= 1.0384; 

U=exp(-(X.-2+Y.~2)) ; U=U/max(max(abs (U) ) ) *A; 
for nn=l :max_iteration 
Uold=U; 

L00U=ifft2(-(KX.-2+KY.~2) . *f f t2 (U) )+3* (cos (X) . ~2+cos (Y) . ~2) . *U + U."3; 
MinvU=ifft2(fft2(U) . /(c+KX . ~2+KY. ~2) ) ; 
mu=sum(sum(L00U. *MinvU) )/sum(sum(U. *MinvU) ) ; 
U=U+if f t2 (f f t2 (L00U-mu*U) . / (KX . ~2+KY . ~2+c) ) *DT ; 
U=U/max (max ( abs (U) ) ) * A ; 

Uerror(nn)=sqrt(sum(sum(abs(U-Uold) . ~2))*dx*dy) ; Uerror(nn) 
if Uerror(nn) < error_tolerance 
break 

end 

end 



This code as well as Matlab codes for other examples can be downloaded at 
www . cems . uvm . edu/ ~ j yang / P ublication . ht m 



References 

[1] N.K. Efremidis, S. Sears, D.N. Christodoulides, J.W. Fleischer and M. Segev, "Discrete soli- 
tons in photorefractive optically induced photonic lattices," Phys. Rev. E 66, 046602 (2002). 

[2] Y.S. Kivshar and G.P. Agrawal, Optical solitons: from fibers to photonic crystals, Academic 
Press, San Diego, 2003. 



21 



[3] J. Yang and Z. Musslimani, "Fundamental and vortex solitons in a two-dimensional optical 
lattice," Opt. Lett. 28, 2094 (2003). 

[4] Th. Busch and J.R. Anglin, "Mossbauer effect for dark solitons in Bose-Einstein condensates," 
http://arxiv.org, Preprint cond-mat 9809408. 

[5] M.L. Chiofalo, S.Succi, and M.P. Tosi, "Ground state of trapped interacting Bose-Einstein 
condensates by an explicit imaginary-time algorithm," Phys. Rev. E 62, 7438 (2000). 

[6] L.D. Carr and Y. Castin, "Dynamics of matter-wave bright soliton in an expulsive potential," 
Phys. Rev. A 66, 063602 (2002). 

[7] H.B. Keller, "Numerical solution of bifurcation and nonlinear eigenvalue problems," in P. H. 
Rabinowitz, ed., Applications of Bifurcation Theory, pp. 359-384; Academic Press, 

[8] J. P. Boyd, "Why Newton's Method is Hard for Travelling Waves: Small Denominators, KAM 
Theory, Arnold's Linear Fourier Problem, Non-Uniqueness, Constraints and Erratic Failure," 
Math. Comput. Simul. 74, 7281 (2007). 

[9] J. Yang, "Internal oscillations and instability characteristics of (2+1) dimensional solitons in 
a saturable nonlinear medium," Phys. Rev. E. 66, 026601 (2002). 

[10] P.G. Kevrekidis, K.O. Rasmussen, and A.R. Bishop, "Localized excitations and their thresh- 
olds," Phys. Rev. E 61, 4652 (2000). 

[11] V. I. Petviashvili, "Equation of an extraordinary soliton," Plasma Physics, 2, 469, (1976). 

[12] D.E. Pelinovsky and Yu. A. Stepanyants, "Convergence of Petviashvili's iteration method 
for numerical approximation of stationary solutions of nonlinear wave equations," SIAM J. 
Numer. Anal. 42, 1110 (2004). 

[13] M.J. Ablowitz, Z.H. Musslimani, "Spectral renormalization method for computing self- 
localized solutions to nonlinear systems", Opt. Lett. 30, 21402142 (2005). 

[14] T.I. Lakoba and J. Yang, "A generalized Petviashvili iteration method for scalar and vector 
Hamiltonian equations with arbitrary form of nonlinearity" , J. Comp. Phys. 226, 1668-1692 
(2007). 

[15] J.J. Garcia-Ripoll and V.M. Perez-Garcia, "Optimizing Schrodinger functional using Sobolev 
gradients: Applications to quantum mechanics and nonlinear optics," SIAM J. Sci. Comput. 
23, 1316 (2001). 

[16] W. Bao and Q. Du, "Computing the ground state solution of Bose-Einstein condensates by 
a normalized gradient flow," SIAM J. Sci. Comput. 25, 1674 (2004). 

[17] V.S. Shchesnovich and S.B. Cavalcanti, "Rayleigh functional for nonlinear systems," available 
at http://www.arXiv.org, Preprint nlin. PS/0411033. 

[18] J. Yang and T.I. Lakoba, "Universally-convergent squared-operator iteration methods for 
solitary waves in general nonlinear wave equations", Stud. Appl. Math. 118, 153-197 (2007). 



22 



[19] J. Douglas and H.H. Rachford, "On the numerical solution of heat conduction problems in 
two and three space variables," Trans. Amer. Math. Soc. 82, 421 (1956). 

[20] S.E. Koonin, Computational Physics, Sec. 7.4; Addison- Wesley, Redwood City, 1986. 

[21] N. G. Vakhitov and A. A. Kolokolov, "Stationary solutions of the wave equation in the 
medium with nonlinearity saturation," Izv. Vyssh. Uchebn. Zaved. Radiofiz. 16, 1020 (1973) 
[Radiophys. Quantum Electron. 16, 783 (1973)]. 

[22] C. K. R. T. Jones, "Instability of standing waves for nonlinear Schrodinger-type equations, 
Ergod. Theory Dynam. Syst. 8, 119 (1988). 

[23] M. Grillakis, "Linearized instability for nonlinear Schrodinger and Klein-Gordon equations," 
Comm. Pure Appl. Math. 41, 747 (1988). 

[24] M. Weinstein, "Lyapunov stability of ground states of nonlienar dispersive evolution equa- 
tions," Comm. Pure Appl. Math. 39, 51 (1986). 

[25] W.A. Strauss, M. Grillakis, J. Shatah, W. Strauss, "Stability theory of solitary waves in the 
presence of symmetry. I", J. Funct. Anal. 74, 160C197 (1987). 

[26] W.A. Strauss, "Existence of solitary waves in higher dimensions," Comm. Math. Phys. 55, 
149 (1977). 

[27] H. Berystycki and P.L. Lions, "Nonlinear scalar field equations I — Existence of a ground 
state," Arch. Rat. Mech. Anal. 82, 313 (1983). 

[28] J. Yang and Z. Chen, "Defect solitons in photonic lattices," Phys. Rev. E. 73, 026609 (2006). 

[29] R. Horn and C. Johnson, Matrix Analysis, Cambridge University Press, New York, 1991. 

[30] M. Struwe, Variational Methods: Applications to Nonlinear Partial Differential Equations 
and Hamiltonian Systems, 3rd ed., Springer, 2000. [Specifically, see the paragraph just below 
Theorem B.4 on page 246.] 

[31] W.A. Strauss, Nonlinear Wave Equations, CBMS Reg. Conf. Ser. Math, vol.73, AMS, Prov- 
idence, RI, 1989. 

[32] E.L. Ince, Ordinary differential equations, Dover, New York, 1956. [Specifically, see the 
Theorem stated at the end of Sec. 10.61.] 

[33] J. L. Bona, P. E. Souganidis, and W. A. Strauss, "Stability and instability of solitary waves 
of Korteweg-de Vries type," Proc. Roy. Soc. London Ser. A 411, 395 (1987). 

[34] R. L. Pego and M. I. Weinstein, "Eigenvalues, and instabilities of solitary waves," Phil. Trans. 
Roy. Soc. London Ser. A 340, 47 (1992). 

[35] E.A. Ostrovskaya and Y.S. Kivshar, "Photonic crystals for matter waves: Bose-Einstein 
condensates in optical lattices," Optics Express 12, 19 (2004). 

[36] N.K. Efremidis, J. Hudock, D.N. Christodoulides, JW. Fleischer, O. Cohen, and M. Segev, 
"Two-Dimensional Optical Lattice Solitons," Phys. Rev. Lett. 91, 213906 (2003). 



23 



[37] Z. Musslimani and J. Yang, "Self-trapping of light in a two-dimensional periodic structure," 
J. Opt. Soc. Am. B. 21, 973 (2004). 

[38] S.V. Manakov, V.E. Zakharov, L.A. Bordag, A.R. Its, and V.B. Matveev, "Two-dimensional 
solitons of the Kadomtsev-Petviashvili equation and their interaction," Phys. Lett. A 63, 205 
(1977). 

[39] X.P. Wang, M.J. Ablowitz, and H. Segur, "Wave collapse and instability of solitary waves of 
a generalized Kadomtsev-Petviashvili equation," Physica D 78, 241 (1994) 



24 



(a) (b) 



u 




-10 -5 5 10 -10 -5 5 10 

X X 



Figure 1: (a) A solitary wave with nodes in the focusing nonlinear Schrodinger equation (40) 
= 1.2689, P = 1), for which the AITEM (16)-(18) diverges for any At; (b) A nodeless solitary 
wave in the defocusing nonlinear Schrodinger equation (41) (// = 3.5069, P = 1), for which the 
AITEM (16)-(18) converges if the stepsize condition At < At max in Theorem 1 is met. 
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Figure 2: Solitary waves with nodes in the defocusing nonlinear Schrodinger equation (42) and 
their L\ spectra, (a) The power diagram; (b) two solitary waves with powers P = 3 and 10; (c, d) 
spectra of L\ for these two waves. The AITEM (16)-(18) converges for the one with P = 10 when 
the stepsize restriction in Theorem 1 is met. 
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Figure 4: (a) A solitary wave in the two-dimensional NLS equation (66) with amplitude one; 
error diagrams of the AITEM (A.N.) and the Petviashvili method for this solitary wave. 




Figure 5: (a) A solitary wave in Eq. (67) with Vq = 3 and P = 3. (b) Error diagrams of 
AITEM and AITEM (A.N.) [both with c = 0.7, At = 1] for this solitary wave. 
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Figure 6: (a) A solitary wave in the KP equation (68) with \i = 1; (b) error diagrams of the AITEM 
(with c = 1.4, At = 1.7), AITEM (A.N.) (with c = 1.4, At = 1.4) and the Petviashvili method for 
this solitary wave. 
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